home *** CD-ROM | disk | FTP | other *** search
/ Shareware Grab Bag / Shareware Grab Bag.iso / 090 / bode.arc / BODE.PAS
Pascal/Delphi Source File  |  1985-09-23  |  8KB  |  213 lines
  1. PROGRAM BODE ;
  2.  
  3. { Program to calculate the frequency response plot (bode diagram)
  4.   for a transfer function provided to the program in terms of
  5.   numerator and denominator polynominal coefficients of the
  6.   Laplace transfer function. The order of the polynominals must
  7.   be less than or equal to 20.
  8.  
  9.   First revision completed by Mike Secord on 10/23/84
  10.          The program was written to try to optimize the execution speed.
  11.   Second revision -- fix bugs          11/15/84
  12.   Revision 2.1    -- display program description when executed  8/34/85 }
  13. CONST
  14.       VERSION='VERSION  2.1' ;
  15.       MAXSTEPS=100;
  16.       NATLOG=2.302585093;
  17.  
  18. TYPE
  19.      NUMFREQ=ARRAY[1..MAXSTEPS] OF REAL;
  20.      COEF_ARRAY = ARRAY [0..20] OF REAL ;
  21.      NAME = ARRAY [1..11] OF CHAR ;
  22.  
  23. VAR  FREQVAL:NUMFREQ;
  24.      NSTEP,I:INTEGER;
  25.      M : INTEGER ;    { ORDER OF NUMERATOR POLYNOMIAL }
  26.      N : INTEGER ;    { ORDER OF DENOMINATOR POLYNOMIAL}
  27.      PN : COEF_ARRAY ; { ARRAY OF NUMERATOR POLYNOMIAL COEFFICIENTS }
  28.      PD : COEF_ARRAY ; { ARRAY OF DENOMINATOR POLYNOMIAL COEFFICIENTS }
  29.      INTYPE : NAME ; { NAME FOR TYPE OF INPUT (DENOM OR NUM) }
  30.      KGAIN, W, MAGN, MAGD, PHASEN, PHASED, MAGNITUDE, LOGMAG, PHASE : REAL ;
  31.      ANS : CHAR ;
  32.  
  33. PROCEDURE TRANFUNC(INTYPE :NAME ; VAR O : INTEGER; VAR C : COEF_ARRAY) ;
  34.  {procedure to input the transfer function (pole and zero) information
  35.  for the BODE program.}
  36.  
  37. VAR
  38.    L : INTEGER ;
  39.  
  40. BEGIN
  41.    WRITELN ;
  42.    WRITE(' Enter the order of the ',INTYPE,' polynomial (20 Max) : ') ;
  43.    READLN(O) ;
  44.    WRITELN ;
  45.    WRITELN(' Enter the coefficients of the ',INTYPE,' polynomial in ') ;
  46.    WRITELN(' desending order (PLEASE ENTER ZERO COEFFICIENTS ALSO).  ') ;
  47.    WRITELN(' ',INTYPE,'  coefficients : ') ;
  48.    FOR L := O DOWNTO 0 DO
  49.    BEGIN
  50.       WRITE('  S',L,': ') ;
  51.       READ(C[L]) ;
  52.       IF( (O-L) MOD 4 = 0) AND (L <> O) THEN WRITELN ;
  53.    END ;
  54.    WRITELN ;  WRITELN ;
  55. END ;
  56.  
  57. PROCEDURE CALCFREQ(VAR FREQ:NUMFREQ;VAR V:INTEGER);
  58. {calculate series of frequencies to be used for bode plot}
  59. VAR     INITIALFREQ:REAL;
  60.         FINALFREQ:REAL;
  61.         LOGDECADE_STEP:INTEGER;
  62.         STARTFREQ:REAL;
  63.         ENDFREQ:REAL;
  64.         TRUNCFREQ:INTEGER;
  65.         DIFF:REAL;
  66.         LOGSTEP:REAL;
  67.         LOGFREQ:REAL;
  68.  
  69. BEGIN {input initial frequency,final frequency,and number of steps per decade}
  70.         WRITELN ;
  71.         WRITE(' Type initial frequency desired in output: ');
  72.         READLN(INITIALFREQ);
  73.         INITIALFREQ := INITIALFREQ - INITIALFREQ/1000.0 ; {Allow for trunc}
  74.         WRITE(' Type final frequency desired in output: ');
  75.         READLN(FINALFREQ);
  76.         FINALFREQ := FINALFREQ + FINALFREQ/1000.0 ; {Allow for truncation}
  77.         WRITE(' Type no. of log steps desired for each decade of output: ');
  78.         READLN(LOGDECADE_STEP);
  79.         {calculate log of input frequencies and frequency steps}
  80.         STARTFREQ:=(1/NATLOG)*LN(INITIALFREQ);
  81.         ENDFREQ:=(1/NATLOG)*LN(FINALFREQ);
  82.         TRUNCFREQ:=TRUNC(STARTFREQ);
  83.         DIFF:=STARTFREQ-TRUNCFREQ;
  84.         LOGSTEP:=1/LOGDECADE_STEP;
  85.         LOGFREQ:=((TRUNC(DIFF/LOGSTEP))/LOGDECADE_STEP)+TRUNCFREQ;
  86.         V:=0;
  87.         {calculate series of frequencies and put into array}
  88.         WHILE LOGFREQ<=ENDFREQ DO
  89.             BEGIN
  90.                 V:=V+1;
  91.                 FREQ[V]:=EXP(NATLOG*LOGFREQ);
  92.                 LOGFREQ:=LOGFREQ+LOGSTEP;
  93.             END;
  94. END;
  95.  
  96. PROCEDURE CALC(ORDER : INTEGER; C : COEF_ARRAY; VAR MAG, ANGLE : REAL;
  97.                W : REAL) ;
  98. { Procedure to calculate the magnitude and phase response for the numerator
  99.   or denominator as called. The total response will be calculated in the
  100.   main body of the program. }
  101.  
  102.    FUNCTION CMULT(W,K : REAL; EXPON : INTEGER) : REAL ;
  103.    { This function calculates the w terms to the respective powers as called
  104.    and was faster than using EXP(U*LN(A)).}
  105.  
  106.    VAR
  107.       P : INTEGER ;
  108.       TEMP2, MPLICND : REAL ;
  109.  
  110.    BEGIN
  111.        MPLICND := W ;  TEMP2 := W ;
  112.        CASE EXPON OF
  113.           0     : CMULT := K ;
  114.           1     : CMULT := W*K ;
  115.           2..20 : BEGIN     { 20 needs to be changed for more than 20th order}
  116.                      FOR P := 2 TO EXPON DO TEMP2 := TEMP2*MPLICND ;
  117.                      CMULT := K*TEMP2 ;
  118.                   END ;
  119.        END ; {case}
  120.    END ; {function CMULT }
  121.  
  122. VAR
  123.    TEMP, RCMULT, ICMULT : REAL ;
  124.    I1 : INTEGER ;
  125.  
  126. BEGIN
  127.    RCMULT := 0.0 ;  ICMULT := 0.0 ;
  128.    FOR I1 := 0 TO ORDER DO
  129.    BEGIN
  130. { NOTE: These sets for the case statement have to be added to
  131.         in increments of 4 if the procedure is
  132.         to calculate the response of a polynominal > 20. }
  133.      CASE I1 OF
  134.        0,4,8,12,16,20 : RCMULT := RCMULT + CMULT(W,C[I1],I1) ;
  135.        1,5,9,13,17    : ICMULT := ICMULT + CMULT(W,C[I1],I1) ;
  136.        2,6,10,14,18   : RCMULT := RCMULT - CMULT(W,C[I1],I1) ;
  137.        3,7,11,15,19   : ICMULT := ICMULT - CMULT(W,C[I1],I1) ;
  138.      END ; {case}
  139.    END ; {for}
  140.    MAG := SQRT(RCMULT*RCMULT + ICMULT*ICMULT) ;
  141.    IF (RCMULT = 0.0) AND (ICMULT > 0.0)  THEN  ANGLE := PI/2.0 ;
  142.    IF (RCMULT = 0.0) AND (ICMULT < 0.0)  THEN  ANGLE := -PI/2.0 ;
  143.    IF (ICMULT = 0.0) AND (RCMULT < 0.0)  THEN  ANGLE := PI ;
  144.    IF RCMULT > 0.0 THEN  ANGLE := ARCTAN(ICMULT/RCMULT) ;
  145.    IF (RCMULT < 0.0) AND (ICMULT > 0.0) THEN
  146.       ANGLE := ARCTAN(ICMULT/RCMULT) + PI ;
  147.    IF (RCMULT < 0.0) AND (ICMULT < 0.0) THEN
  148.       ANGLE := ARCTAN(ICMULT/RCMULT) - PI ;
  149. END ; { procedure ccalc }
  150.  
  151. {*************** MAIN PROGRAM ****************************************}
  152.  
  153. BEGIN
  154.    CLRSCR ;
  155.    GOTOXY(15,8) ;
  156.    WRITE('********** PROGRAM BODE **********') ;
  157.    GOTOXY(26,9) ;
  158.    WRITE(VERSION) ;
  159.    GOTOXY(1,11) ;
  160.    WRITELN(' Program to calculate the frequency response plot (Bode diagram)');
  161.    WRITELN(' for a transfer function provided to the program in terms of');
  162.    WRITELN(' numerator and denominator polynominal coefficients of the');
  163.    WRITELN(' Laplace transfer function. The order of the polynominals must');
  164.    WRITELN(' be less than or equal to 20.');
  165.    WRITELN ;
  166.    WRITELN(' Enter data when prompted (Do NOT forget the zero coefficients)');
  167.    WRITELN ;  WRITELN ;
  168.    INTYPE := ' numerator ' ;
  169.    TRANFUNC(INTYPE,M,PN) ;
  170.    INTYPE :='denominator' ;
  171.    TRANFUNC(INTYPE,N,PD) ;
  172.    WRITELN ;
  173.    WRITE(' Enter the gain coefficient (Enter 1 if there is none) : ') ;
  174.    READLN(KGAIN) ;
  175.    CALCFREQ(FREQVAL,NSTEP);
  176.    WRITELN ;  WRITELN ;
  177.    WRITELN('  FREQUENCY     MAGNITUDE     MAGNITUDE DB     PHASE DEG ') ;
  178.    WRITELN ;
  179.    FOR I:= 1 TO NSTEP DO
  180.    BEGIN
  181.       W := FREQVAL[I]*2.0*PI ;
  182.       CALC(M,PN,MAGN,PHASEN,W) ; {Calculate the numerator respsonse}
  183.       CALC(N,PD,MAGD,PHASED,W) ; {Calculate the denominator response}
  184.       MAGNITUDE := KGAIN* MAGN/MAGD ;
  185.       LOGMAG := (20.0/NATLOG)*LN(MAGNITUDE) ; {Magnitude in DB}
  186.       PHASE :=(PHASEN - PHASED)*180.0/PI ; {Phase in degrees}
  187.       WRITELN(' ',FREQVAL[I]:12,'   ',MAGNITUDE:12,'   ',LOGMAG:12,
  188.               '   ',PHASE:12);
  189.    END;
  190.  
  191.    {Part of program to print to printer on request}
  192.    WRITELN;
  193.    WRITE(' Would you like the output directed to the printer (Y or N) : ') ;
  194.    READLN(ANS) ;
  195.    IF ANS ='Y' THEN
  196.    BEGIN
  197.       WRITELN(LST) ;  WRITELN(LST) ;
  198.       WRITELN(LST,'  FREQUENCY     MAGNITUDE     MAGNITUDE DB     PHASE DEG ');
  199.       WRITELN(LST) ;
  200.       FOR I:= 1 TO NSTEP DO
  201.       BEGIN
  202.          W := FREQVAL[I]*2.0*PI ;
  203.          CALC(M,PN,MAGN,PHASEN,W) ; {Calculate the numerator respsonse}
  204.          CALC(N,PD,MAGD,PHASED,W) ; {Calculate the denominator response}
  205.          MAGNITUDE := KGAIN* MAGN/MAGD ;
  206.          LOGMAG := (20.0/NATLOG)*LN(MAGNITUDE) ; {Magnitude in DB}
  207.          PHASE :=(PHASEN - PHASED)*180.0/PI ; {Phase in degrees}
  208.          WRITELN(LST,FREQVAL[I]:12,'   ',MAGNITUDE:12,'   ',LOGMAG:12,
  209.               '   ',PHASE:12);
  210.       END; {FOR}
  211.    END ; {END PRINTER IF}
  212. END. {BODE}
  213.